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We use the finite-entanglement scaling of iMPS to explore supposedly infinite order transitions. 
This universal method may have lower computational costs than finite-size scaling. To this end we 
study possible MPS based algorithms to find the ground states of the transverse ANNNI model in a 
spin chain with first and second neighbour interactions and frustration. The ground state has four 
distinct phases with transitions of second order and one of supposedly infinite order, Kosterlitz- 
Thouless transition. To explore phase transitions in the model, we study general quantities like 
the correlation length, entanglement entropy and the second derivative of the energy with respect 
to the external field, and test the finite-entanglement scaling. We propose a scaling ansatz for the 
correlation length of a non critical system in order to explore infinite order transitions. This method 

provides considerable less computational costs compared to the finite-size scaling method in Ref. 

gl and quantities obtained by applying fixed boundary conditions (like domain wall energy in Ref. 

31) are omitted. The results show good agreement with previous studies of finite-size scaling using 
DMRG. 



I. INTRODUCTION 

Recently, numerical simulations of many-body quantum systems based on matrix product states 
(MPS) and their generalizations have called increasing attention as they offer an efficient way to 
study properties of their ground states and even thermal states. One of the main reasons of their 
popularity is their applicability to-in principle-any Hamiltonians describing only local interactions 
on a spin chain and their possible generalizations to higher dimensions. The use of the so-called 
projected entangled pair states (PEPS) to approximate ground states of such kind of Hamiltonians 
defined in higher dimensions (e.g. on a square lattice) is considered very promising. (See e.g. Ref.Q 
and references therein for a review of the topic.) These methods can provide efficient simulation 
of systems when other methods such as quantum Monte Carlo break down, e.g. in case of certain 
frustrated systems. It means for example that one can circumvent the so-called sign problem present 
in the Monte Carlo simulation. (Note that there exist systems for which there is no other method 
available, not considering exact diagonalization which works only for small system sizes. One such 
system is e.g. the frustrated XX model. Ref. [l^ ) 

Obviously one of the main goals of investigating ground states is to explore phase transitions and 
characterize them as precisely as possible. Despite the fact that by definition, the MPS can only 
describe exponentially decaying correlations, one can obtain accurate results for certain quantities 
of phase transitions like in the Ising model as it has been demonstrated in nearly all of the papers 
presenting an MPS based algorithm. (See e.g. Ref. where it has also been used for an Ising model 
defined on a binary tree.) Models with not only nearest neighbor (NN) interactions and higher than 
second order transitions have not been studied extensively by means of MPS. 

This paper has two goals: Find a suitable efficient algorithm based on MPS to simulate Hamiltonians 
with second neighbor interactions and frustration, as well as demonstrate the use of universal tools, 
especially finite-entanglement scaling to study different types of phase transitions. Its validity arises 
from general facts about critical points and the definition of the MPS, hence it can be applied quite 
universally. As opposed to the usual finite-size scaling it may also have a better scaling for the 
computational costs if using the same maximal bond dimensions for both cases (see below). Here 
the 1-dimensional version of the transverse axial next nearest neighbor Ising (ANNNI) model has 
been chosen which is non-integrable and possesses a quite complex phase diagram, despite being the 
simplest model with NNN interactions. This model is notable due to several facts. As earlier studies 
have shown, unlike the Ising model, it has a critical region, i.e. a 2D parameter region where it 
is critical. A probably infinite order phase transition also takes place (besides second order phase 
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transitions). As mentioned earlier, the detection of an infinite (e.g. Kosterlitz-Thouless) transition is 
a notoriously hard computational challenge, that's why a computationally beneficial method is needed. 
The same holds for exploring the physics of a critical region. 

In order to study such a model one can choose a finite chain of length L with periodic (PBC) or 
open boundary conditions (OBC) and then look at the scaling of some physical quantities depending 
on L and the bond dimension D. It is known that the overall computational costs scales 0{D^) in 
the best case. The latter happens for example when the correlations are very small ^ <C i as it is 
discussed in Ref. [9.], however when ^ « L e.g. it can be as large as 0{D^) and the precise behaviour 
of the scaling as a function of D and L is unclear. 

Note that in this paper, IMPS and imaginary time evolution are used, hence one works directly 
on the infinite lattice and the only scaling parameter is D. In this case the computational cost 
in every small time step At scales as 0{D^). Roughly speaking, convergence can be reached if: 
exp(iAi?) = exp(nAtAi?) ^ 1 with the IS.E energy gap and n the number of time steps. Approaching 
to a critical point, vanishes as a function of the effective correlation length ^d, and at the same 
time we also have a scaling relation like oc (see below). With this in mind, it can be estimated 
how the overall computational time scales as a function of D and one sees why it is slower than 0{D^) 
near to a critical point. For instance, if /S.E oc oc D^"*^, a > then n should scale as 0{D°'^) 
that makes the overall scaling 0{D^'^'^^). Previous work about finitc-D scaling can be found e.g. in 
and tia]. 

Now the idea is to study the scaling of physical quantities not as a function of the system size, 
but simply as a function of the bond dimension D. Instead of finite-size scaling this can be called 
as 'finite-D' or 'finite-entanglement' scaling as a MPS with a given D can contain a finite amount of 
entanglement. As it will be argued in Section IV. the overall computation cost of finite-D scaling can 
be considerably reduced compared to that of the usual finite-size scaling. In fact, this turns out to be 
the case for the presented model as well as for the Ising and Heisenberg models, and essentially this 
depends on the central charge. The key point to this computational cost reduction is our off critical 
scaling ansatz Q for the correlation length as a function of D which turns out to be computationally 
more beneficial than that of used in Ref. Q for the entanglement entropy. 

Clearly, our methods can be applied straightforwardly to any other ID model with NN and NNN 
interactions, such as the Ji — J2 Heisenberg model, but can also be adapted for the 2D version using 
PEPS. Despite the polynomial computational time scaling, unfortunately its exponent is so high in 
the 2D case that in practice one can apply only very small bond dimensions providing only a few and 
ambiguous results that are not yet sufficient for numerical analysis. 

II. THE ID TRANSVERSE ANNNI MODEL 

Consider the following Hamiltonian on an infinite chain of spin-I/2 particles 

H = - ^( JlO-f Cr,f_^i + J20'iCr^^2 + f^'^i) — Hising + HmNN- 
i 

This can be considered as the 1 dimensional version of the lattice model with axial next nearest 
neighbor interactions (ANNNI). If the next nearest neighbor (NNN) interaction term describes an- 
tiferromagnetic interactions (J2 < 0) and Ji > then the system is frustrated. Note that basically 
two parameters determines the physical properties of model e.g. k = J2/J1 and h/Ji . Changing the 
basis of every second spin appropriately by rotating around the axis x leads to Z^* — Z^', so one 
can switch from an antiferromagnetic NN interaction to a ferromagnetic one and vice-versa, leaving 
the NNN interaction unchanged. However, there is no rotation switching all couplings to ferromag- 
netic, thus the sign of the NNN interactions is crucial and besides the magnetic field of h/\Ji\ the 
ratio K = J2/\Ji\ plays the key role in determining the physical properties. From now on let us set 
Ji = I > and J2 <0 yielding competing ferro-and antiferromagnetic NN and NNN couplings. 
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III. MPS STATES AND ALGORITHMS 

The MPS representation of an infinite, translational invariant chain of d-dimensional systems looks 
hke 

d 

|*)= Tr(ArA^^.)kl,S2...) 

Si,S2.. = l 

with site independent sets of matrices {A^js of dimension D x D. The normahzation constraint 
requires that the largest eigenvalue of the transfer matrix E — has to be unique and 

equals to 1. If it's degenerate (as it happens e.g. for the state 4* = co ttt ■•■ + ci Hi ■■■) then it can 
cause troubles in some types of the algorithms that work with unique leading eigenvectors (Ref. 0, 
Ref. Q). Nevertheless, in order to overcome this one can apply tricks like appropriate rotations of 
some basis vectors, but due to the lots of computations like leading eigenvectors, inverses and squares 
of matrices these algorithms turn out to be costly. Another representation Ref. [3] uses the Schmidt 
coefficients in the diagonal matrices A' explicitly as follows 

I*) = ■••r(*)'"'A(')r('+i)'^'+iA('+i)...|s,)|s,+i)... 

with the normalization constraints 

Tr(A«)^ ^ 1 

J2 (rW-^-)^(Aw)'rW-^- = Id, J2 (rW'^')^(A(*-^))'r«'^' = 

s=l s=l 

in the general case when the dimensions of A(^-i) and A(*) can be different. The infinite time evolving 
block decimation (iTEBD) Ref. [4] uses this description and in case of NN interactions two neighboring 
sites are updated after the action of e^''*^*'»+i. One builds up a larger matrix occupying the MPS 
matrices of the two sites i and i + 1 and decomposing them through a singular value decomposition 
to obtain the new matrices. A normalization step can be interpreted as a zero time update: we just 
build a larger matrix and decompose it. The details can be found in Ref. [2]- 

Several ways exist to do the updates in the case at hand. One option is the use of 'superspins' 
( made of 2 neighboring spins). To this end let us join together every second neighboring pairs of 
spins to form block spins with 2x2 = 4 dimensional basis as ti+i) |ii+2, ^i+s) = |si)|si+2) with 
ti G {0,1} and Si G {0,1,2,3}. In this way we get a translational invariant Hamiltonian with only 
NN interactions acting between these 4 dimensional block spins so we can use the simplest algorithms 
designed for this case but with d — 4 and a bit more complicated NN term. 

In our scheme, for example the operator erf (E> cf+i £ Af4(C) acting on the sites i, i + 1 in the original 
lattice is interpreted as an operator acting on the single superspin labelled by i, while 1^ (Ei <^i+i ® 
l'+3 ^ -^^16(C) describes an interaction between the superspins at the position i and i + 2. 
Similarly trf (g) l^+i ® crf^j ® lj+3 6 -^16 (C) is a NN interaction between the superspins, however it 
describes a NNN coupling between the original spins. With this in mind one can easily construct the 
operators 

exp(eicr,f(7,f+i + e2cr,f o-f+a) 

in the new basis. 

We make use of the 'superspins' and the update method described in Ref.[l]. The computation 
time in every step scales as 0{D^) but as opposed to the previous method like in Ref.0, one does not 
need to compute neither leading eigenvectors nor inverses or squares of non diagonal matrices hence 
it is faster and there is no problem with degenerate eigenvalues. 
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IV. FINITE-ENTANGLEMENT SCALING TO DETECT PHASE TRANSITIONS 

In order to characterize phase transitions using iMPS one needs to study how some special quantities 
scale as a function of the applied bond dimension D. Let us denote the entropy of entanglement 
computed from the IMPS by S]j and the correlation length ^^i. It is known that in the critical points 
Sd converges logarithmically with corrections of order 0(1/ log as 

Sd = ^ log D (1) 
with c being the central charge and k only depends on c for large enough D as 

12+1 



The scaling of the correlation length looks like 

cx (3) 

again with corrections of 0(1/ log D), see Ref. 1^ and [ll[ for more details. 

In practice however, one should go in the other way round: Testing the scaling relations one wants 
to decide whether the system is critical or not by giving an estimate for the correlation length. The 
closer we are to the critical point the more delicate the problem due to the high correlation length. 
One has to decide if the limits S = Soo = liniD^oo So and ^ = ^oo = limD^cx) are finite or not. {S 
and ^ converge to the exact values of the infinite system as the Trotter errors tend to zero.) It can be 
especially hard for example in the vicinity of a Kosterlitz-Thouless transition when the inverse of the 
true correlation length tends to zero very quickly as we are approaching the critical point. As in 
finite-size scaling let us assume some scaling function for as 

logCzj =fclogi? + .9(ACoo) (4) 

Clearly, g{D, ^oo) ~ log(^oo/Z)*'') for large -D if Ci? ^ Coo < oo. It enables us to define f{x) — 
g{D,^ao) with X = D'^/^oo, for which f{x) — > if x ^» and f{x) ^ — log(a;) for large x. Note that 
besides these limit properties we will not make any further observations about the form of /, just 
use a simple ansatz what will be justified by the numerics to be a good indicator for the non-critical 
behaviour. ( See below the similar ansatz function for the finite-size case). So for instance, the following 
simple ansatz meets these requirements 

f{x)^-\og{x + e-°'^)+ const, a > 0. (5) 

Near the critical regime ^ is big but not infinite and for a not big enough D it behaves very similarly 
to ^ and one cannot distinguish it from the critical scaling. So the natural question arises: For a 
given and large but finite ^ near the critical regime, at least how big D should be applied in order 
to rule out the critical scaling with some reasonable accuracy? Essentially it depends on the ratio of 
logD&ndfif). 

In order to decide whether a given point is critical or not one should test the linearity of log 
against log _D as a null hypothesis. If it can be rejected then one can fit a trial function according to 
(in and ([5]) to obtain an estimate for ^oo- 

Note thatQ resembles of the finite-size scaling of the entropy of the half chain of length L applied 
in Ref. Q 

5(L) = |(logL + /(LC')) (6) 

with the same ansatz for f{x) but now one has to plug x = -Z^/^oo- Note that this simple ansatz ([5]) 
for f{x) has the right limits but for intermediate values it is not correct, however it turns out to be 
good approximation for S{L) and a good indicator for non-criticality by the numerics. One can obtain 
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FIG. 1: Schematic phase diagram 

(universal) series expansions for f{x) in the small and large limit theoretically, see Ref. [T^. One can 
also confront this with our ansatz which incorporates the to limits. 

Of course here, for a given L, the applied bond dimension (or the number of states kept in each 
DMRG iteration) should be "high enough". 

In our case however, we only have the bond dimension as an input parameter which simplifies the 
picture. Besides, the main advantage is as follows: Having the same ^oo for the infinite system, for 
the same accuracy one has considerably smaller computational cost if the scaling exponent fc > 1. 
(This will turn out to be the case in the present model as the central charges are found to be c = 1 
implying k w 1.34 and c = 1/2 implies k w 2.03.) In any case, the time cost of one update step scales 
as 0{LD^) in the finite size case as opposed to 0{D^) for the infinite case. Now, if one takes the 
same maximal bond dimension D for both cases, then for the finite- scaling one needs to also apply 
some smaller D-s which makes the overall time costs scale as 0{D^) (at most, if one passes along 
1, 2...D — \,D). On the other hand in case of finite-size scaling the time costs oc 0{L'^D^). In order to 
detect off critical scaling with the same accuracy one expects the necessary maximal length L in the 
finite-size case must be bigger than the maximal D in the finite-D case, if fc > 1. Nevertheless, this 
can be estimated by studying the scaling anzatz ^ and ©. In terms of D this would result in an 
overall effective scaling higher than 0{D^), in the finite-size case, obtained by the assumption L (x D. 

V. RESULTS 

Figure [Tj shows the phase diagram as it has been found in earlier studies by simulations Ref. Q and 
perturbative analysis Ref.Q. We try to explore the phase boundaries by looking at the correlation 
lengths as a function of the parameters and the second derivative of the energy. A Kosterlitz-Thouless 
type transition is conjectured between the floating and the paramagnetic(PM) phase where the deriva- 
tives of the energy are useless indicators so the expected divergence of the correlation length is checked 
only. In the critical points we also check the finite entanglement scaling laws and find very good agree- 
ment with the theoretical results. Moreover, in some special points we calculate the critical charges 
as well as the correlation scaling exponents. 

We intend to determine some points of the line of phase transitions between the distinct phases and 
check these relations as well as the finite-entanglement scaling and compute the central charge. There 
are three lines of phase transitions, two of them are of second order while an infinite order transition 
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is expected for the third one which needs a special care. We also try to detect the floating phase when 
K < and |k| is large, since its existence for arbitrary large |k| is still an open question. 

A. FM-PM transition 

In case of the Ising model, in the vicinity of the critical point as a function of the magnetic field, 
the correlation length and the magnetization diverges as 

(x \h - hct" , M(x\h-hcf 

and the energy and its first derivative are continuous but its second derivative diverges at h — he- 
Using D = 40 we get 0.535 < he < 0.538 by looking at the extrema of the finite derivative d^E/dh^. 
At the transition corresponding to k = —0.25 we find the same behaviour for the energy as a function 
of h and a very good agreement with v = 1 as in the simple Ising case: v = 1.02 comes from the linear 
fit on the figure below. 
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J2/J^=-.25 
D=32 




□ 

0.6- 



-5.6 -5.4 -5.2 -5 -4.8 -4.6 -4.4 -4.2 -4 -3.8 

log|h-hJ 

FIG. 2: Estimation of the scaling exponent at k = —0.25 paramagnetic-ferromagnetic transition 



B. Antiphase-floating phase transition 

Again a second order phase transition is suggested by the numerics but with a central charge c = 1 
of the floating phase as opposed to the Ising case with c = 1/2. See figure 

C. Floating phase-paramagnetic phase transition 

First, we check the power-law scaling of ^ assumed for critical points as in ([3]). If this power-law 
scaling can be rejected statistically, then for a given h and k, S^oo is computed by the best fit using (g]) 
and finite-D scaling. (See the inset in FIG|H) We observe that fits well using the ansatz function 
with weak dependence of the parameter a which can be set to 1. 

In fact, the numerics indicate a higher order transition here, with continuous derivatives of the 
energy. Probably a Thouless-Kosterlitz type transition takes place as suggested by the plot vs h 
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FIG. 3: Entanglement entropy scaling inside the floating phase at k — —0.75, h — 0.4 up to D = 100. 
Central charge c = 1 is suggested according to ([T]) and ([2}. 




Ql , , , , , , , 1 

0.34 0.36 0.38 0.4 0.42 0.44 0.46 0.48 0.5 

FIG. 4: Floating phase-paramagnetic phase transition illustrated with the inverse correlations plotted for 
three different bond dimensions. The inner plot shows l/log(^oo) vs h. Here, for a fixed h, ^oo is computed 
by fitting (g]). 

since £,{h)~^ oc exp( ^ ) could be fitted ii h > he, a > 0. At the point k = —0.75 around h — 0.45 
very close to the floating/PM transition we need about Z? > 80 to point out the non power-law scaling 
(see. FIG. [S]). Otherwise, where we cannot reject the power-law scaling hypothesis, even though a 
finite value for ^oo obtained by fitting ^ is rejected. In these cases one cannot rule out the critical 
scaling statistically. We have got anyway quite ambiguous results for which -in this case-can be 
regarded as a consequence of the errors of the simulation. As in Ref. [1], in FIG. |3]we plot the inverse 
correlation length along the line k = —0.75 near to the supposed Kosterlitz-Thouless transition and 
find similar results with those calculated by finite DMRG and studying domain wall energies. Having 
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FIG. 5: against D in the two sides of the floating-paramagnetic transition. It can be seen in the log-log 
scale how it bends down from the straight line for higher D, implying a non critical point a,t h = 0.5. On the 
other hand, at h = 0.4 the critical scaling hypotesis ([3} cannot be rejected with the applied maximal bond 
dimension 

our division for the magnetic field, the critical value can be estimated roughly as 0.42 < he < 0.44. 
Note that here, the study of model specific quantities have been avoided. 

D. Large frustration 

At At = —10 up to the the applied maximal bond dimension Dmax = 80 and the division for the 
magnetic field, the numerics do not seem to indicate the presence of the fioating phase because one 
sees a sharp peak plotting ^ against h. See FiglBl 

Note that in Ref. Q the existence of the floating phase up to k = — 5 has been indicated but 
the DMRG method has become imprecise for higher frustrations, however with our new method, the 
floating phase can be suspected even with Dmax — 50 as shown in FIG. \7\ 

VI. SUMMARY AND CONCLUSIONS 

The phase diagram of the ID transverse ANNNI model is studied. We have coupled together 
neighboring sites of the chain in order to get NN interactions between the 4-dimensional superspins 
and applied the simple iTEBD with the new NN Hamiltonian. In practice, the iTEBD algorithm seems 
to be more beneficial than others based on the transfer matrix and its leading eigenvectors. At some 
typical points of the phase diagram we have studied the phase transitions and estimated the critical 
quantities and found good agreement with previous DMRG results coming from finite-size analysis 
and model specific quantities. We have also checked finite-entanglement scaling and its breaking to 
localize a supposedly infinite order transition between the floating and the paramagnetic phase, as 
well as pointed out the logarithmic divergence of the entanglement entropy inside the floating phase. 
With the proposed scaling ansatz for f one can identify the critical regime with less computational 
effort in comparison with flnite-size scaling, hence for the same results lower bond dimensions are 
sufficient. For a high frustration parameter n — —10 the numerics do not indicate the presence of 
the fioating phase with the applied maximal bond dimension D = 80 and resolution of the magnetic 



9 




FIG. 7: The correlation length against the magnetic field using different bond dimensions at k = —5 indicates 
an extended critical interval. 



field Ah ~ 0.08. However, we can suspect the existence of the floating phase for k = —5 already 
with D <= 50. The question of its existence for large frustrations requires a more detailed study 
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with higher bond dimensions, but the apphed methods throughout the paper are quite universal and 
apphcable to other frustrated spin systems as welL 
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